#include "funcLogistic.h"

using namespace QPCRModel;
funcLogistic::funcLogistic(){
    a = 1;
    b = 1;
    x0 = 1;
    y0 = 0;
    d = 1;
}

double funcLogistic::f(double x){
    return y0 + a / (1 + pow(x/x0, -b));
}

double funcLogistic::g(double y){
    return x0 * pow(a/(y-y0) - 1, -1/b);
}

double funcLogistic::loss(double x, double y){
    return y - y0 - a/(1/pow(x/x0, b) + 1);
}

double funcLogistic::loss(double x, double y, const Eigen::VectorXd &X){
    return y - X[3] - X[0]/(1/pow(x/X[2], X[1]) + 1);

}

double funcLogistic::loss2(double x, double y){
    return pow(loss(x, y), 2);
}

double funcLogistic::loss2(double x, double y, const Eigen::VectorXd &X){
    return pow(loss(x, y, X), 2);
}

double funcLogistic::diff1(double x){
    return (a*b)/(x0*pow(1/pow(x/x0, b) + 1, 2)*pow(x/x0, b+1));
}

double funcLogistic::diff2(double x){
    return (2*a*b*b)/(x0*x0*pow(1/pow(x/x0, b) + 1, 3)*pow(x/x0, 2*b + 2)) - (a*b*(b + 1))/(x0*x0*pow(1/pow(x/x0, b) + 1, 2)*pow(x/x0, b + 2));
}

Eigen::VectorXd funcLogistic::jacobian(double x, double y){
    Eigen::VectorXd res(4);
    res <<  -1/(1/pow(x/x0, b) + 1),
            -(a*log(x/x0))/(pow(1/pow(x/x0, b) + 1, 2)*pow(x/x0, b)),
            (a*b*x)/(x0*x0*pow(1/pow(x/x0, b) + 1, 2)*pow(x/x0, b + 1)),
            -1;
    return res;
}

Eigen::VectorXd funcLogistic::jacobian2(double x, double y){
    Eigen::VectorXd res(4);
    // 好像用不上....不算了
    return res;
}

void funcLogistic::setV(const Eigen::VectorXd &v){
    a = v[0];
    b = v[1];
    x0 = v[2];
    y0 = v[3];
}

Eigen::VectorXd funcLogistic::getV(){
    Eigen::VectorXd res(4);
    res << a, b, x0, y0;
    return res;
}

double funcLogistic::getLinearStageFunc(double *lk, double *lb){

    double x_fdm = x0 * pow((b-1)/(b+1), 1/b);
    double y_fdm = this->f(x_fdm);
    *lk = this->diff1(x_fdm);
    *lb = -(*lk) * x_fdm + y_fdm;
    return *lk;
}


double funcLogistic::curvity(double x, double norm_coef){
    if(norm_coef <= 0){
        norm_coef = a;
    }
    double res = ((2*norm_coef*b*b)/(x0*x0*pow(1/pow(x/x0, b) + 1, 3)*pow(x/x0,2*b + 2)) - (norm_coef*b*(b + 1))/(x0*x0*pow(1/pow(x/x0, b) + 1, 2)*pow(x/x0, b + 2)))/pow((norm_coef*norm_coef*b*b)/(x0*x0*pow(1/pow(x/x0, b) + 1, 4)*pow(x/x0,2*b + 2)) + 1, 1.5);
    return res;
}

int funcLogistic::getParamNum(){
    return 4;
}

std::pair<double, double> funcLogistic::getRange(int i){
    if(i == 0){
        return std::pair<double, double>(0, INFINITY);
    }else if(i == 1){
        return std::pair<double, double>(0, 100);
    }else if(i == 2){
        return std::pair<double, double>(0, 60);
    }else if(i == 3){
        return std::pair<double, double>(0, INFINITY);
    }else{
        return std::pair<double, double>(-INFINITY, INFINITY);
    }
}

void funcLogistic::init(double SPE, double EPE, int cycleNum, double ymin, double ymax){
    a = 0;
    b = 0;
    x0 = 0;
    y0 = 0;
    d = 1;

    // 根据findExp2函数，SPE和EPE定义为线性期两侧斜率小于2的点
    if(EPE > cycleNum){
        a = (ymax - ymin) * 2;
    }else{
        a = (ymax - ymin);
    }
    b = 12;
    x0 = (SPE + EPE) / 2;
    y0 = ymin;


}
